Exact relativistic kinetic theory of an electron beam-plasma system: 
hierarchy of the competing modes in the system parameter space 
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Besides being one of the most fundamental basic issues of plasma physics, the stability analysis 
of an electron beam-plasma system is of critical relevance in many areas of physics. Surprisingly, 
decades of extensive investigation had not yet resulted in a realistic unified picture of the multidi- 
mensional unstable spectrum within a fully relativistic and kinetic framework. All attempts made 
so far in this direction were indeed restricted to simplistic distribution functions and/or did not aim 
at a complete mapping of the beam-plasma parameter space. The present paper comprehensively 
tackles this problem by implementing an exact linear model. We show that three kinds of modes 
compete in the linear phase, which can be classified according to the direction of their wavenumber 
with respect to the beam. We then determine their respective domain of preponderance in a three- 
dimensional parameter space. All these results are supported by multidimensional particle-in-cell 
simulations. 



Relativistic electron beams moving through a collision- 
less background plasma are ubiquitous in a variety of 
physical systems pertaining, e.g., to inertial confinement 
fusion py], the solar corona [1], electronic pulsar winds 
j^, accreting black hole systems active galactic nu- 
clei or 7-ray burst sources 0. Regardless of its set- 
ting, the electron beam-plasma system has been known 
to be prone to collective processes since the very early 
studies on the so-called two-stream instability [7] , which 
is of electrostatic nature and therefore characterized by 
wave and electric field vectors both parallel to the beam 
flow direction. Later on, another class of unstable modes 
was discovered: usually referred to as "filamentation" [8|] 
or Weibel modes [9], they are mostly electromagnetic, 
purely growing and develop preferentially in the plane 
normal to the beam. Finally, unstable modes propa- 
gating obliquely to the beam were investigated [lO| and 
found to rule the system in case of cold and diluted rela- 
tivistic electron beams 11, 12, T3|. To date, the few ki- 
netic approaches investigating the whole two-dimensional 
(2-D) unstable spectrum have been restricted to non- 
relativistic energy spreads [14] or to waterbag-like distri- 
bution functions of questionable validity at high energy 
spreads [15]. Among the salient features of the oblique 
modes revealed by the latter studies are their mostly elec- 
trostatic character and usually efficient interaction with 
both the beam and plasma components. Yet, because of 
various limitations, all the models considered so far have 
failed to provide a complete picture of the hierarchy of 
the competing unstable modes as a function of the sys- 
tem parameters. One should stress that addressing this 
long-standing issue, which is the main goal of the present 
work, is a critical prerequisite to understanding the non- 
linear aspects of the beam-plasma interaction. Among 
these are the instability-induced collective stopping and 
scattering, whose level is expected to depend closely on 
the nature (i.e., electrostatic or electromagnetic) of the 
underlying processes [11, [13, H E, 0, [HI, [H, [11, [H . 



model implemented here involves, for both the counter- 
streaming beam and plasma populations (hereafter iden- 
tified by the subscripts h and p), unperturbed distribu- 
tion functions in the form of drifting Maxwell- Jiittner 
functions SSEzl 



exp [-/ic.(7(p)-/^aP2y)] ' (1) 



In contrast to past studies, the fully relativistic kinetic 



which allow for arbitrary energy spreads and drifts. Here 
a = (b^p) stands for the beam or plasma component. 
Pa = {Py/l) is the normalized ^/-aligned mean drift ve- 
locity, 7c, the corresponding relativistic factor and /j^a = 
mc^/kBTa the normalized inverse temperature of each 
electron component. All momenta are normalized by 
mc. K2 denotes a modified Bessel function of the sec- 
ond kind. Current neutralization is assumed, that is, 
ni)l3i) + riplSp = 0, where and Up are the beam and 
plasma mean densities, respectively. From now on, the 
ions form a fixed neutralizing background and collisions 
are neglected. It is noteworthy that the above model dis- 
tribution function is provided a thermodynamically con- 
sistent derivation from first principles in Ref. [H, [2^, 
where it is shown to maximize the specific entropy for 
fixed values of each species' total momentum and energy. 
In addition to its thermodynamic grounds, this distribu- 
tion function is helpful in gaining insight into the stability 
properties of smooth relativistic distributions, as opposed 
to the commonly used waterbag distributions 
severely flawed by the neglect of Landau damping. The 
evolution of the initially homogeneous and unmagnetized 
system is governed by the relativistic Vlasov-Maxwell set 
of equations. Following a routine procedure, the general 
dispersion equation for any orientation of the wave vec- 
tor is derived [13, HH and solved numerically. Details 
of the numerical procedure will be presented elsewhere. 
Introducing the total electron plasma frequency cjg, and 
normalizing the complex frequency a; and wave vector 
k by uoe and uOe/c respectively, there remain four inde- 
pendent variables: the beam mean relativistic factor 75, 
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FIG. 1: (color online). Normalized growth rate in the {kx,ky) 
plane for n^/rip = 1, 75 = 1.2, Tb = 500 keV and Tp = 5 keV. 
Isocontours are linearly spaced from 0.01 to 0.07. 

the beam and plasma temperatures T5 and Tp, and the 
density ratio rib /rip. A typical map of the growth rate 
S = ^uj is displayed in Fig. [T] for n^/n^ = 1, 75 = 1.2, T5 
= 500 keV and Tp = 5 keV. The three aforementioned in- 
stability classes are clearly visible for this configuration. 
For wave vectors aligned with the beam, the two-stream 
instability peaks for ky ^ 0.5 whereas, in the perpendic- 
ular direction, the fastest growing filamentation mode is 
found for kx ~ 0.5. Overall, though, the 2-D unstable 
spectrum is here dominated by an oblique mode located 
at (kx.ky) - (0.5,0.5). 

The three instability classes do not share the same sen- 
sitivity to the beam temperature, drift and density. This 
yields a non-trivial dependence of the growth rates on 
the system parameters, hence a varying hierarchy be- 
tween the instability classes. Assuming from now on a 
fixed plasma temperature Tp = 5 keV, it is possible to de- 
termine the regions of predominance of each instability 
class in the {nb/rip^jb^Tb) space. The surfaces that de- 
limit regions governed by different instability classes are 
displayed in Fig. and colored according to the local 
maximum growth rate. Points located between the plane 
75 = 1 and the left surface define systems dominated 
by the two-stream instability, while those located be- 
tween the right surface and the plane rib/ rip = 1 pertain 
to filamentation-ruled systems. Oblique modes prove to 
govern the rest of the parameter space. 

The two-stream instability is seen to prevail in the 
whole non-relativistic range of the beam drift energy 
(75 — 1 <C 1), as well as in weakly relativistic systems 
with hot enough beams. The reason for the latter fea- 
ture is as follows: first, filamentation is strongly weak- 
ened for weakly energetic beams because its growth rate 
is proportional to the beam velocity. Being thus left with 
the two-stream and oblique modes, we note that the for- 
mer are handicapped by the relativistic increase in the 
longitudinal inertia, while the latter are more sensitive 
to the beam temperature. As a consequence of these ef- 
fects, a system characterized by rib /rip = 0.1, 75 = 1.5 
and Tb = 500 keV turns out to be mostly subject to 
the two-stream instability, as indicated in Fig. [2](a). We 



have checked this prediction by means of a 2D particle- 
in cell (PIC) simulation using the massively parallel code 
CALDER [32]. As expected, the plasma density profile, 
displayed in Fig. Efb) at the time when the instabil- 
ity starts to saturate, exhibits modulations preferentially 
along the drift direction. By repeating the Fourier anal- 
ysis performed in Ref. [l5|, we have verified that the 
simulated k-resolved growth rates closely agree with lin- 
ear theory. 

Let us now consider the filamentation-to-oblique tran- 
sition, which is found to take place for dense and rel- 
ativistic beams. The shape of the corresponding fron- 
tier is here more involved than the previous one as 
it stems from a balance between the three system pa- 
rameters. The filamentation growth rate increases with 
the beam density and decreases with temperature more 
rapidly than the oblique mode growth rate. For a cold 
system [13], this growth rate scales like /^b/y^, which 
reaches a maximum for 75 = a/3. Still in the cold 
limit, the filamentation-to-oblique frontier reaches a min- 
imum rib /rip ^ 0-53 for 75 ~ 3. The resulting bound- 
ary proves to be mostly determined by 75 for dense, 
cold and weakly relativistic beams. In the relativistic 
and ultra-relativistic regimes, the main parameter is the 
beam density, although high enough beam temperatures 
always end up favoring oblique modes. This result goes 
against the conventional belief that relativistic systems 
with rib /rip = 1 are governed by filamentation 0, [3] • 
We show here that this behaviour holds true only as long 
as the beam is not too hot. This is illustrated in Fig. 
[2jc,d) which compares the simulated plasma density pro- 
files for two values of the beam temperature, respectively 
below and above the filamentation-to-oblique boundary. 
The other parameters are rib /rip = 1 and 75 = 1.5. In 
accord with linear theory, filamentation modes are seen 
to prevail when T5 = 100 keV, whereas oblique modes 
clearly take over when T5 = 2 MeV, yielding a character- 
istic knitted-like pattern. 

When crossing one of the transition surfaces displayed 
in Fig. [2](a), a phase velocity discontinuity may occur 
as regards the dominant mode. During the oblique/two- 
stream transition, the ky component of the most unsta- 
ble mode remains almost constant while the kx compo- 
nent steadily decreases down to zero [Fig. [31(a), curve 
2]. The oblique/two-stream transition is therefore con- 
tinuous, notably as regards the resulting phase veloc- 
ity change. During the oblique/filament at ion transition 
[Fig. [3](b)], by contrast, the kx component hardly varies, 
whereas the ky component evolves as follows. Let us start 
from a system ruled by an oblique mode with ky ^ O and 
7^ [Fig. [3](b), curve 1]. As the system moves to 
the filamentation regime, the ky component of the domi- 
nant mode remains almost constant while its growth rate 
gradually weakens. When the latter eventually becomes 
smaller than that of the main filamentation mode (with 
ky = 3?cj = 0), the system experiences a sudden change 
regarding its dominant wave phase velocity, which then 
switches from a finite value to zero. 

So as to gain insight into the field generation within 
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FIG. 2: (color online). Left: Hierarchy of the unstable modes in the (n^/np, 7^, T^) parameter space for Tp—b keV. The left 
surface delimits the two-stream-dominated domain (at low 75) and the oblique-mode-dominated domain, whereas the right 
surface delimits the filamentation-dominated domain (at high nh/rip) and the oblique- mode-dominated domain. 
Right: Plasma density profiles at the end of the linear phase as predicted by 2D PIC simulations run with three different sets 
of parameters: nh/rip = 0.1, 7^ = 1.5 and Th = 500 keV (b); nh/rip = 1, 7^ = 1.5 and Th = 100 keV (c); nh/rip = 1, 75 = 1.5 
and = 2 MeV (d). In all cases, Tp = 5 keV. In agreement with linear theory, the three resulting patterns evidence regimes 
dominated by two- stream, filament at ion and oblique modes, respectively. 
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FIG. 3: (color online). Typical evolution of the fastest grow- 
ing mode (red arrow) for the oblique/two- stream (a) and 
oblique/filamentation (b) transitions, of which only the for- 
mer is continuous. 



a more realistic geometry, we have performed a three- 
dimensional PIC simulation of a diluted beam-plasma 
system defined by 77,5 /n^ = 0.1, 75 = 3, T5 = 50 keV 
and Tp = 5 keV (Fig. gj). In agreement with Fig. [2l^a), 
the system is initially governed by oblique modes, as con- 
firmed by the beam and plasma profiles at uoet = 80. 
Later on {uOet = 160), though, the system gets ruled by 
two-stream modes. A rough, quasilinear-like argument 
supporting this transition can be given assuming that 
the beam distribution function has retained a Maxwell- 
Jiittner form. At uOet = 160 the best fit is obtained for 
76 ^ 1.6 and T5 ~ 200 keV. Because of the low beam den- 
sity, we can reasonably neglect the changes in the plasma 
distribution function. Figure [2](a) then indicates that the 
system is indeed dominated by two- stream modes. After 
saturation of all electrostatic modes, the system eventu- 
ally reaches a filamentation regime, owing to a remaining, 



high enough anisotropy, as already observed in Ref. [lH • 
Yet, this transition cannot be explained in light of the 
present linear model because the distribution functions 
then strongly depart from the form ([T]). 

Let us conclude by making two important remarks. 
Firstly, given our model distribution function ([T]), com- 
plete linear stabilization of the system can never be 
achieved in the whole parameter space. In particular, 
the maximum filamentation growth rate can be shown 
to scale like in the large T5 limit. Secondly, the 

oblique/filamentation boundary in the plane T5 = 1 keV 

— 1/3 

behaves like rib/up ~ 1 — 0.867^ when 75 1. As 
a result, the filament ation-ruled region is all the more 
squeezed against the plane Ub/up = 1 when the beam en- 
ergy is high. This also implies that the ultra-relativistic 
regime is asymptotically oblique unless rib /rip = 1- This 
could bear important consequences in a number of astro- 
physical scenarios which involve relativistic factors up to 
10^ and beyond [35l, 
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Oblique Two-Stream Filamentation 

FIG. 4: (color online). 3D PIC simulation initialized with nh/up — 0.1, 7b = 3, = 50 keV and Tp = 5 keV: isosurfaces of the 
beam and plasma density profiles at coet = 80, 160 and 560. The system runs through three successive phases, each governed 
by a distinct instability class. 
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